temp = load('temp.dat');

Temp = temp(:, 1);
Time = temp(:, 2);

co = cos(2*pi*Time/72);
si = sin(2*pi*Time/72);

[seasonb, seasonbint, seasonres] = regress(Temp, [ones(size(co)), co, si]);

prefilt = [[seasonres(3:end); 0; 0], 2*[seasonres(2:end); 0], 3*seasonres, 2*[0; seasonres(1:(end-1))], [0; 0; seasonres(1:(end-2))]]/9;
postfilt = sum(prefilt, 2);

%Empirical transfer
%te = abs(fft(postfilt)./fft(seasonres)).^2;
%Theoretical transfer
tt = abs(fft([1/9, 2/9, 3/9, 2/9, 1/9, zeros(1, 347)])).^2;

subplot(3, 1, 1)
plot(seasonres, '-k', 'LineWidth', 2)

set(gca, 'Box', 'off', 'FontSize', 20, ...
         'XLim', [0, 360], 'YLim', [-0.6, 0.6], ...
         'XTick', 0:72:288, 'YTick', -0.6:0.3:0.6, ...
         'XTickLabel', 0:4, 'TickDir', 'out')
xlabel('Time (days)', 'FontSize', 20)


subplot(3, 1, 2)
plot((1:176)/352, tt(1:176), '-k', 'LineWidth', 2)

set(gca, 'FontSize', 20, 'XLim', [0, 0.5], 'XTick', 0:0.1:0.5, 'YLim', [0, 1])
xlabel('Frequency', 'FontSize', 20)

subplot(3, 1, 3)
plot(postfilt, '-k', 'LineWidth', 2)

set(gca, 'Box', 'off', 'FontSize', 20, ...
         'XLim', [0, 360], 'YLim', [-0.6, 0.6], ...
         'XTick', 0:72:288, 'YTick', -0.6:0.3:0.6, ...
         'XTickLabel', 0:4, 'TickDir', 'out')
xlabel('Time (days)', 'FontSize', 20)

set(gcf, 'Position', [200, 100, 1400, 900])